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We study the jamming transition in a model of elastic particles under shear at zero temperature. 
The key quantity is the relaxation time r which is obtained by stopping the shearing and letting 
energy and pressure decay to zero. At many different densities and initial shear rates we do several 
such relaxations to determine the average r. We establish that r diverges with the same exponent 
as the viscosity and determine another exponent from the relation between r and the coordination 
number. Though most of the simulations are done for the model with dissipation due to the motion 
of particles relative to an affinely shearing substrate (the RDq model), we also examine the CDq 
model, where the dissipation is instead due to velocity differences of disks in contact, and confirm 
that the above-mentioned exponent is the same for these two models. We also consider finite size 
effects on both r and the coordination number. 

PACS numbers: 63.50.Lm, 45.70.-n 83.10.Rs 


I. INTRODUCTION 

Granular materials, supercooled liquids, and foams are 
examples of systems that may undergo a transition from 
a liquid-like to an amorphous solid state as some con¬ 
trol parameter is varied. It has been hypothesised that 
the transitions in these strikingly different systems are 
controlled by the same mechanism [ 1 ] and the term jam¬ 
ming has been coined for this transition. Much work 
on jamming has focused on a particularly simple model, 
consisting of frictionless spherical particles with repul¬ 
sive contact interactions at zero temperature [2]. The 
packing fraction (density) (j) of particles is then the key 
control parameter. Many investigations have focused on 
jamming upon compression, and jamming by relaxation 
from initially random states [2-4]. Another, physically 
realizable and important case, is jamming upon shear 
deformation. This has been modeled with elastic parti¬ 
cles both with a finite constant shear strain rate 7 [5-11], 
and by quasistatic shearing [4, 12, 13], in which the sys¬ 
tem is allowed to relax to its local energy minimum after 
each finite small strain increment. A nice method to do 
shearing simulations of hard disks has also recently been 
developed[14]. 

Several open questions remain in spite of much stud¬ 
ies of the jamming models under steady shear. Central 
among them is an understanding of the mechanisms be¬ 
hind jamming, a question that has been addressed, for 
the case of hard disks, in several papers by Wyart and 
co-workers[14-16]. A related question is what details of 
the models that are important for the universality class. 
It has earlier been claimed[ll] that a more realistic model 
for the dissipation—where the dissipation is due to the 
velocity differences between disks in contact, the CDg 
model—gives a different critical behavior than the sim¬ 
pler RDq model in which the dissipation is against an 
affinely shearing substrate. Evidence agains this claim 
has recently been given in[17], but much work remains 
to clarify other aspects of the various models that are 
relevant for different physical systems close to jamming. 


In this work we perform large scale simulations to 
determine the relaxation time—a quantity whose diver¬ 
gence, we will argue, lies behind the jamming transition. 
We do that by first shearing at a steady shear rate and 
then stopping the shearing and letting energy and pres¬ 
sure decay to zero; the relaxation time is the time con¬ 
stant of this exponential decay. We also determine a re¬ 
lated time—the dissipation time—which is the time scale 
of the initial decay just after stopping the shearing. We 
characterize the dependencies of these relaxation times 
on both distance from (below) jamming and the initial 
shear rate. We then motivate a direct relation between 
the relaxation time and the lowest vibrational frequency 
of Lerner et al.[14]. Following Lerner et al.[I4] we deter¬ 
mine the contact number z in the absence of rattlers. We 
then find that the relaxation time depends algebraically 
on the distance to the isostatic contact number, and de¬ 
termine the exponent for this divergence. Most of our 
simulations are for the simpler RDq model (see below) 
but we also do the same kind of analysis for the CD 
model, and confirm[17] that these two models appear to 
behave the same. We then turn to two effects that are 
related to the finite system sizes: We first show that the 
ordinary arithmetic averaging can sometimes give un¬ 
expected effects, and then examine how the number of 
particles in the simulations affects the spread in contact 
number and relaxation time. 

The organization of this paper is as follows: In Sec. II 
we describe our numerical methods and give a brief sum¬ 
mary of some earlier results that are used throughout the 
paper. In Sec. Ill we first introduce our two key quan¬ 
tities and discuss their differences and similarities. We 
then discuss the relation to the vibrational frequencies 
in a model of hard disks[I4]. Also following Ref. [14], 
we demonstrate a direct relation to the contact number 
and show that the determined exponent is the same for 
CDq as for RDq. We also consider the finite size effects. 
In Sec. IV we finally discuss our results, relate them to 
earlier works, and make some comments. Sec. V gives a 
short summary. 
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II. MODEL AND SIMULATIONS 

A. Simulations 

Following O’Hern et al. [2] we use a simple model of 
bi-disperse frictionless soft disks in two dimensions with 
equal numbers of disks with two different radii in the 
ratio 1.4. Length is measured in units of the diameter 
of the small particles, dg- With the distance between 
the centers of two particles and dij the sum of their radii, 
the interaction between overlapping particles is V (rij) = 
{ej2)S'^j with the relative overlap 6ij = 1 — nj/dij. We 
use Lees-Edwards boundary conditions [18] to introduce 
a time-dependent shear strain 7 = tj. With periodic 
boundary conditions on the coordinates Xi and yt in an 
L X L system, the position of particle z in a box with 
strain 7 is defined as = {xi +jyi,yi)- The simulations 
are performed at zero temperature. 

We consider two different models for energy dissipa¬ 
tion. The CD model (CD for “contact dissipation”) is 
the model introduced by Durian for bubble dynamics in 
foams [19], and was also used by Tighe et al. [11]. Here 
dissipation occurs due to velocity differences of disks in 
contact, 

fcD,* = -kd - Vj), V, = r,. (1) 
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In the second model, RD—“reservoir dissipation”—the 
dissipation is with respect to the average shear flow of a 
background reservoir, 

fRD.i = -kd{^y^ - VR(ri)), VR(ri) = yz/iX. (2) 

RD was also introduced by Durian [19] as a “mean-field” 
[20] approximation to CD, and is the model used in many 
early works on criticality in shear driven jamming [5, 14, 
20, 21]. In both cases the equation of motion is 

= + (3) 

We are here interested in the overdamped limit, mz —>■ 0 
[19]. In the RD model it is straightforward to perform 
simulations with m = 0. In the CD model we take m = 1 
which, for the shear rates we are using, turns out to be 
small enough to be in the overdamped limit. We take 
e = 1 and fcd = 1. The unit of time is tq = dskd/e. 

We focus most of our effort, using longer simulation 
runs at lower shear rates, for the model RDq, but we also 
give results for the model CD for comparison. We use 
N = 65536 particles, and shear rates down to 7 = 10“® 
and integrate the equations of motion with the Heuns 
method with time step At = 0.2ro. 

B. Background 

The present paper focuses on the behavior of the 
above-mentioned models just below (pj. We here sum¬ 
marize a few results that are important in the following. 


The jamming transition is a zero-temperature transi¬ 
tion from a liquid to a disordered solid upon the increase 
of density. An excellent way to probe this transition is 
to look at the resistance to shearing. Since the defining 
property of a liquid is that it is a material that cannot 
sustain a shearing force, a finite shear stress, cr, in the 
limit 7 —>■ 0 is a clear sign of a solid phase. Within the 
liquid, i.e. aX (p < pj, the approach to jamming is seen 
in the rapid increase of the viscosity, 7 = < 7 / 7 ; numerical 
evidence suggest that it diverges algebraically, 

r](p, 7 0) = a/j ~ {pj - p)~^. (4) 

Another quantity that clearly signals the transition is 
the pressure and the pressure equivalent of the viscosity, 
dp = p/ii which similarly diverges with the exponent j3, 

r]p{p,p ^ 0) =p/j ^ {pj - p)~/^. (5) 

Since p ^ 6 whereas the interaction energy is A ~ <5^, the 
energy diverges with the exponent 2/3, 

\im E/j'^ ^ {pj - p)~‘^^. (6) 

7^0 

Equations (4) and (5) for cr and p, should hold very close 
to pj, but since the dimensionless friction, p = ajp = 
dldp^ has a strong ^dependence, Eqs. (4), (5) clearly 
give only the leading divergence of 7 and rjp, and are 
not exact expressions that hold over any finite density 
interval. To handle this one needs to include corrections 
to scaling[22, 23] by writing 

0/jr^{Pj-P)-^[l + coiPj-Pri ( 7 ) 

for the observables cr and p. 

In simulations of soft particles the data will depend 
on the shear rate, 7 , which may be considered a relevant 
scaling variable. This suggests a scaling assumption as 
in critical phenomena[5]. With 5p = p — pj, 

0{Sp,P) = b-y/'^goiSpb^/'',W), ( 8 ) 

where b is typically considered to be a length rescaling 
factor, though it can be chosen arbitrarily. With b = 
\5p\~'', specializing to 5p < 0, the scaling relation for 
O/p becomes 

0{5p,P)/P=\5p\-^^''-y^go{il\5pr)- ( 9 ) 

In the 7 —>■ 0 limit go{x —>■ 0) = const, together with 
Eq. (7) leads to the identification pi = zv — y. 

Corrections to scaling are included by generalizing 
Eq. ( 8 ) to 

O = b-y/''[go{5pb^/\ib^) + b-‘^ho{5pb^/\ib^)], (10) 

An analysis based on this kind of approach[22] gave 
P = 2.77(20) whereas a related approach in terms of an 
effective density [24] gave the very similar /3 = 2.58(10). 
Other recent values in the literature from simulations 
are /3 = 2 . 2 [ 21 ], and a recent theoretical work gives 
P = 2.77[16]. 
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III. RESULTS 
A. Measured quantities 

1. The relaxation time 

One of the hallmarks of the jamming transition is a 
diverging time scale. It has been common to measure 
this time scale implicitly by measurement of a diverging 
transport coefficient like 77 or rjp. In this section, how¬ 
ever, we measure such a time scale by looking directly at 
the relaxation of the system from an initial shear driven 
steady state to the zero-energy state obtained after the 
shearing is turned off. We thus make use of a two-stage 
process: In the first stage the system is driven at steady 
shear with a constant shear rate 7 , in the second stage 
the shearing is stopped but the dynamics is continued 
which makes the system relax down to a minimum en¬ 
ergy. As the simulations discussed here are at densities 
somewhat below (j)j, the final state is always a state of 
zero energy, and after a short transient time, energy and 
pressure decay exponentially to zero. The relaxation time 
for a single relaxation is denoted by ti , 

p{t) - exp(-t/Ti). 

A few such relaxations at different densities are shown 
in Fig. 1. In each case the relaxation time is determined 
from the data with p{t) < 10 “^, where the decay is ex¬ 
ponential to an excellent approximation. As we will see 
below the relaxation time depends on the shear rate ap¬ 
plied before the relaxation and we will let 7 )—which 

thus depends on both </> and 7 —denote the average re¬ 
laxation time from about 10-100 such relaxations. 

Fig. 2(a) which is 7 ) versus (j) for several differ¬ 
ent shear rates, clearly suggests that r diverges at the 
jamming transition. The figure also illustrates the shear 
rate dependence; r gets bigger for larger 7 which means 
that the system driven at higher shear rates needs longer 
time for reaching the zero-energy state. The reason for 
this behavior is maybe not entirely obvious, but one can 
at least say that the opposite behavior—that the decay 
were faster for a higher initial shear rate—would be very 
counterintuitive. Recall that this is the shear rate before 
the relaxation step; the relaxation itself is performed with 
7 = 0. 

Note also that this relaxation time is a different quan¬ 
tity from the quantity with the same name in the context 
of supercooled liquids. In supercooled liquids the parti¬ 
cles’ motion is due to the non-zero temperature, whereas 
the motion in the present context is due to the relaxation 
of the potential energy. 

2. Dissipation time 

As a complement to the relaxation time, which is de¬ 
termined from the final decay of the pressure, we also 



FIG. 1. Examples of the pressure relaxation at different (p. 
The figure shows the pressure relaxation after the shearing has 
been switched off. The preceeding shearings were performed 
at very low shear rates in order to stay close to the linear 
region; the densities and the initial shear rates were {(j>, 7 ) = 
(0.8340,10“®), (0.8380,10"®), (0.8400, 5 x 10"'^), (0.8408, 2 x 
10“®), (0.8416, 10“®). To determine the relaxation times, r, 
we fit pressure to an exponential decay, only using data with 
p < 10“T 


introduce the “dissipation time” Tdissj which is defined 
from the initial decay rate, just after the shearing has 
been turned off. For this quantity there is however no 
need to study the actual relaxations; at any moment the 
relaxation rate for the energy may be determined from 
the energy together with the dissipating power, giving 
t' = E/P^iss- In steady shear we may equate the dis¬ 
sipated power with the input power Pin = U(T 7 , which 
gives t' = Ej{a'j) for the average dissipation time. As we 
want a quantity that may be directly compared to r—i.e. 
the decay time for pressure rather than the decay time 
for energy—we note that p S whereas E 6'^ which 

means that p{t) ~ 6 “*/’’ implies E{t) ^ and 

that the two relaxation times differ by a factor of two. 
Our final expression for the dissipation time is therefore 

E , , 

Tdiss=2 —. (II) 

aj 

Figure 2(b) shows Tdiss against (p for several different 
shear rates. Just as for r this quantity also appears to di¬ 
verge as (j) ^ (pj. The 7 -dependence is however different; 
Tdiss decreases with increasing 7 , which means that the 
relative decrease of the energy is bigger in simulations 
at higher shear rates. The different behaviors of r and 
Tdiss is presumably because Tdiss picks up contributions 
from all kinds of decay modes, and the faster modes are 
more excited when the system is driven with a higher 
shear rate. In contrast, t only gets contributions from 
the slowest decay mode. 

Figure 2(c) shows a comparison of t and Tdiss- To elim¬ 
inate effects due to the finite shear rate we only include 
the data at the smallest shear rate and exclude the data 
at p = 0.8420, 7 = 10“® which is clearly away from the 
7—^0 limit. 
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FIG. 2. Relaxation time and dissipation time vs density. 
Panel (a) shows r vs 0 at several different shear rates. The 
data increases rapidly with increasing (f) suggestive of a di¬ 
vergence at (j)j. There is also a clear shear rate dependence, 
T decreases when 7 is decreased towards the hard disk limit, 
7 —>■ 0. Panel (b) which shows raiss vs <j) also increases rapidly 
with 4>j. The shear rate dependence is however the opposite; 
Tdiaa increases with decreasing 7 . Panel (c) shows a compari¬ 
son of r and raiaa which only includes the data with the lowest 
7 (i.e. closest to the hard disk limit), r and raiaa behave essen¬ 
tially the same across this density interval, they are very close 
at the highest density close to <j)J, but the (relative) difference 
increases with decreasing (j>. 




FIG. 3. Divergence of r and raiaa. We here fix (f>j = 0.8433 
and determine /3 by fitting the few points of r and raiaa, re¬ 
spectively, with (j)j — (j> < 0.006 and sufficiently small 7 to be 
close to the hard disk limit. (The points for (f) = 0.8420 and 
7 = 10 ~® appear to be too far from the hard disk limit and 
are not included in the fits.) 


nent /3: 


_ 

"^diss — / . 

a/-f 


[(jij - 4>)~^ 


{(jjj-fl)) 


( 12 ) 


That T diverges in the same way follows from the very 
similar behaviors in Fig. 2(c) but in Sec. IIIB2 we will 
also argue for a direct connection between rjp and r by 
other means. 

Figs. 3 show the determination of P from r and Tdiss- 
The determinations are based on the data points from in 
Fig. 2(c) very close to (j)j, (pj — cj) < 0.006. With only 
a few points with limited precision in a narrow interval 
of (p, it is difficult to do a fit with both /3 and pj as 
free parameters. We therefore instead determine /3 after 
fixing the jamming density to pj = 0.8433[12, 22, 24]. 
The actual fits of r and ruiss are shown in Figs. 3 and give 
similar values for the exponent: /3 = 2.71 and /3 = 2.78 
in good agreement with earlier estimates[16, 22, 24]. 


3. Divergence B. Relations to hard disk simulations 


We are now ready to demonstrate one of the key re¬ 
sults of the present paper, which is that both r and ruiss 
diverge with the exponent p. From the definition of the 
dissipation time in Eq. (11) together with Eqs. (6) and 
(4), it follows directly that ruiss diverges with the expo- 


In this Section we will relate the relaxation time r to 
results from the study of the vibrational modes of sheared 
hard disks [14]. Relations between these two approaches 
are expected since soft disk simulations at sufficiently low 
shear rates give vanishingly small overlaps and therefore 
should behave just as hard disks. 
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1. Relaxation time and the vibrational frequency 


To motivate the relation between the relaxation time 
and the vibrational frequency we consider small displace¬ 
ments Ui from a zero-energy state. Written in terms of 
the vector u, with 2N elements, and the stiffness matrix 
M, such that the force (also a vector with 2N compo¬ 
nents) becomes eAfu, the equation of motion for inertial 
dynamics may be written 




(13) 


(We here consider a finite mass although our work is con¬ 
cerned with the overdamped limit of m —> 0 , only to be 
able to relate to other approaches.) With eigenvalues 
and eigenvectors the force due to a general dis¬ 
placement field, u = J2k becomes eJ2k 

and the ansatz u(t) = sinijjfet gives = 

— {e/m)X^^\ However, below (pj where the number of 
contacts is below the isostatic value there are modes with 
zero energy and iu;fc = 0 , which complicates the analysis. 
From the formalism for shearing of hard disks Lerner et 
al. [14] derived a matrix with the same eigenvalues as Ai 
except for these zero-energy modes. For that matrix the 
lowest frequency, Wmin) is always finite. 

The relaxation may similarly be analyzed in terms of 
small displacements and for overdamped dynamics the 
equation of motion becomes 


kd— = eMu. (14) 

at 

The ansatz of an exponential decay, u(t) = 
exp(-t/rfe) then gives = -{e/kd)X^^^. 

Taken together, Eqs. (13) and (14) give the desired re¬ 
lation between the relaxation time and the vibrational 
frequencies. 


kd -2 

Tk = —. 

m 


(15) 


Our largest —the same as our relaxation time, r—then 
corresponds to the lowest frequency, comin- 

^ - <?„■ ( 16 ) 


Our observation that there is only a single relaxation 
time that controls the decay corresponds well with the 
finding[14] that the lowest frequency in the vibrational 
analysis is an isolated mode. If that were not the case, 
one would expect several decay modes with similar relax¬ 
ation times and that would be seen through a curvature 
in the data in Fig. 1. 


2. Relation to pressure 
The formalism of Ref. [14] gives the relation 

^min Vpi 



0.830 0.835 0.840 


0 


FIG. 4. Comparison of r and r/p which, up to a constant 
prefactor, are expected to behave the same as 7 0. Panel 

(a) shows the raw data, r, and Aprjp, with the constant Ap = 
36. The data for large shear rates (solid triangles) are clearly 
different, but the respective points approach one another as 
7 —>■ 0. The dashed line is /r(</>) ~ (0j — 4>)~^ with /3 = 2.60. 
Panel (b) show the same data, but now divided by frif))- The 
figure clearly suggests that the data should agree in the 7 —>■ 0 
limit. 


to be valid in the hard disk limit. Together with Eq. (16) 
this leads us to expect that r and ijp should behave the 
same in the hard disk limit and Eigs. 4 shows comparisons 
of T and rjp from our soft disk simulations with different 
shear rates. The data clearly approach one another as 
7 —>■ 0 . 


Panel (a) shows t together with Apirjp (where the con¬ 
stant is Ap = 36) against (f) for different 7 . Both quan¬ 
tities do indeed appear to approach the same curve in 
the 7 —0 limit, given by the dashed line, /t(((>) 

{(j)j — Panel (b) which shows the same data, but 

now relative to fT{4')i serves as a strong confirmation of 
the expected equality and gives ample support for the 
expected direct proportionality between t and rip in the 
hard disk limit. Recall that rjp and r are very different 
quantities as the first is determined at constant shear¬ 
ing whereas the second is from the relaxation rate of the 
pressure. 
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C. Contact number 

1. Relaxation time and contact number 

A key result from the study of static packings is that 
jamming in frictionless systems occurs when the coordi¬ 
nation number is z = Ziso = 21 ), which is the number 
needed for mechanical stability. [25] This is however ex¬ 
act only in the absence of rattlers—particles that are not 
locked up at a fixed position as they have less than three 
contacts. To eliminate rattlers we follow Ref. [14] and 
repeatedly remove all particles with less than three con¬ 
tacts. After removing the rattlers, Zi is obtained as the 
average number of contacts of the remaining particles. 

Following Lerner et al. [14] we show the individual de¬ 
terminations, Ti against Szi = Ziso — zi in Fig. 5(a). The 
figure gives strong evidence for an algebraic relation. For 
the vanishing of Szi we introduce Uz, 

5z ^ . (17) 

Together with r ^ {(pj — this gives a relation be¬ 
tween the individual data points ti and zi, 

n ~ (<5zi)-^/“% (18) 

and a fit of our data gives the exponent P/uz = 2.69. 
Since there is a curvature in the data that sets in around 
5zi = 0 . 1 , only data with 5zi < 0.08 were used in the 
fit. This result appears to be especially robust since it is 
obtained from a very simple fit of the raw data with no 
adjustable parameter. (Compare Fig. 3 where a deter¬ 
mination of /3 depends on the correct value of cpj.) Note 
also that there is no need to restrict the data to small 
shear rates of the initial simulation stage. As shown in 
Fig. 5(b) data for different 7 do indeed fall on (or spread 
around) the same line. The explanation for this seems 
to be that both ti and zi are determined from configu¬ 
rations with almost vanishing overlaps, essentially in the 
hard disk limit, independent of the initial shear rate. To¬ 
gether with /3 = 2.70[22] this suggests Uz = 1 whereas 
the somewhat smaller (3 = 2.58 [24] which would imply 
Uz ~ 0.96, means that we cannot exclude the possibility 
that Uz takes on a non-integral value. 

Our result P/uz = 2.69 is in good agreement with 
Ref. [14] who found P/uz = 1/0.38 = 2.63. A more 
recent paper by the same authors[16], however, suggests 
P/uz = 1/0.3 ~ 3.3 (their Fig. 5(c)). This new and lower 
exponent (0.3 < 0.38) is due to a curvature in their data, 
bending over from a larger slope for i5z > 0.1 to this lower 
slope for 5z < 0.1. This bending over at 5zi ~ 0.1 is sim¬ 
ilar to our Fig. 5(a), though the slopes are different. We 
cannot offer any explanation for this difference. (The ef¬ 
fect in Fig. 8 (a) below, which also leads to a larger value 
ol P/uz, doesn’t seem to be applicable in that case.) 

As mentioned above the contact numbers were deter¬ 
mined from the relaxed configurations with almost van¬ 
ishing particle overlaps. To check if it would be possible 
to do a similar analysis of the configurations before the 




FIG. 5. Corresponding values of n and 5zi. Panel (a) shows 
2719 corresponding values of ri and 5z\. Each point is from a 
relaxation that gives both a relaxation time n and a Hnal con¬ 
figuration from which the contact number zi is determined. 
The relaxation time clearly depends algebraically on 5zi —the 
distance to isostaticity. A fit of all data with Szi < 0.08 (1625 
points) gives the exponent Pjuz = 2.69. Panel (b) is a zoom- 
in with a more restricted set of data: p = 0.8412 and four 
different shear rates. This shows that the points for different 
initial shear rates fall on a single curve. 

relaxations, we have also determined the corresponding 
starting values, and to see how the relaxation pro¬ 

cess changes the contact number Fig. 6 shows the final 
contact number, zi against the corresponding starting 
values, zf"*. These data are obtained for p = 0.8412, 
closely below pj, and four different shear rates. From 
the figure we may draw a few different conclusions: (1) 
The contact number always decreases in the relaxation 
process. (2) This change is bigger for larger initial shear 
rates. (3) The final zi decreases slowly with decreas¬ 
ing initial shear rate. (4) The contact number of the 
starting configurations is sometimes above isostaticity, 
_^start ^ whereas zi is always below. This last point 
makes clear that the analyses above, where the approach 
to jamming is seen by zi —>■ Ziso can not be used with 
zstart; obtained from the relaxed configura¬ 

tions that approaches Zjso as jamming is approached. 

2. Analysis of the CDq model 

We have also applied the methods discussed above to 
the CDq model. These results are from a rather lim¬ 
ited number of relaxations and no data very close to 
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FIG. 6 . Change in contact number in the relaxation process. 
The figure shows contact numbers before and after the relax¬ 
ation. The solid line is zi = The configurations are at 

density (f) = 0.8412; the starting configurations are generated 
with four different initial shear rates. Both the initial 
and zi, obtained after the relaxation, are calculated after re¬ 
peatedly removing all particles with less than three contacts. 



FIG. 7. Determination of P/uz for the CDq model. By fitting 
data for Szi < 0.08 to Eq. (18) we determine P/uz = 2.63. 
We note that this is very close to P/uz = 2.69 of the RDq 
model. 

jamming, but they nevertheless give convincing results. 
Fig. 7 shows Ti vs 5zi just as in Fig. 5. The solid line, 
from fitting the data with 5zi < 0.08, gives the expo¬ 
nent P/uz = 2.63. We note that this is very close to 
P/uz = 2.69 of the RDq model which gives support to 
the recent claim[17] that these two models have the same 
critical behavior. To facilitate a direct comparison, the 
fitting line in Fig. 5 is included as a dashed line in Fig. 7. 
The only difference appears to be that the the relaxation 
time for the CDq model is about a factor 1.5 larger than 
for the RDq model, for the same value of 5zi. 


3. Effect of large fluctuations 

Figure 5 above displayed the individual data points 
{ti,Szi), with different symbols for different simulation 
parameters p, 7 . An obvious way to show the same thing 
in a less crowded figure, would be to determine the arith¬ 
metic means of ri and zi for the different sets {p, 7 ). We 
introduce the notation Tq and (Sz)a for these arithmetic 




FIG. 8 . Mean values of ri and Szi determined in two different 
ways. Panel (a) shows the ordinary arithmetic mean values. 
For small Sz these points deviate clearly from the expected al¬ 
gebraic behavior. This phenomenon is due to the large spread 
of the data which appears close to jamming as is also described 
in conjuction with Figs. 9. Panel (b) which shows the geomet¬ 
ric means, Tg, and (Sz)g of the points (ri, Szi) in Fig. 5(a) for 
the same p and 7 . These points obey an algebraic behavior 
with the exponent P/uz = 2.68 in very good agreement with 
the analysis of the individual data points in Fig. 5(a). 

means, (tq is thus just the ordinary average, t.) This 
kind of data is shown in Fig. 8 (a), and it then turns out 
that the averaged data don’t behave quite the same as 
the individual points; the few points at the smallest {Sz)a 
are now clearly off the solid line. The reason for this is 
that the ti for a certain combination of p^ 7 are spread 
over a finite range of 6z and since there is a power law 
relation between r and Sz, if one does the arithmetic av¬ 
erage of this fixed phi data, one gets a point that does 
not lie on the same curve. 

However, it turns out that things work differently— 
all the data fall on the line—when one instead plots the 
geometric means, 

TgiP,!) = exp , (19) 

{Sz)g{p,p) = exp . ( 20 ) 

This data is shown in Fig. 8 (b). 

To illustrate what happens when one averages data 
with a power law relation, Figs. 9 show the behavior of 
arithmetic and geometric means for some points on the 
line y = x~^, on logarithmic and linear scales, respec¬ 
tively. The points labelled “arithmetic” and “geometric” 
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FIG. 9. Illustration of the arithmetic mean and the geometric 
means for some points on the curve y = x~^ . From the figure 
with linear scale in panel (b) it is clear that one cannot expect 
the arithmetic mean to lie on top of the curve. As discussed 
in the text this effect only becomes important in cases where 
the relative variance is sizeable. 

are the respective averages of the open circles in the fig¬ 
ures. In the left panel, which shows the data on loga¬ 
rithmic scales, the arithmetic average is again, just as in 
Fig. 8 (a), clearly off the line. Though this could seem 
surprising, a plot with linear scales as in panel (b) di¬ 
rectly shows that the arithmetic average cannot lie on 
that line. 

This effect is directly related to the big spread in the 
data around the average together with a power differ¬ 
ent from one. With points yi = ya(l + 5i) where ya is 
the arithmetic mean and 5i the relative deviation from 
this mean, the variance is cr^ = (j/^) — (j/)^ = (^?)- 

To second order in the deviations, the geometric mean 
becomes 

Vg = exp((ln[j/a(l + ^»)])) 



FIG. 10. Finite size and the spread of the points 
for (f) = 0.838 and initial shear rate 7 = 10“^. Panel (a) 
which is for three different system sizes shows that 

the points spread considerably more for smaller N. Panel (b) 
shows different quantitative measures of the spread of these 
data. The open circles are s(«i)—the standard deviation of 
zi. Open squares are s(Ti)/ra. (The normalization by Ta is 
to get quantities of the same order of magnitude). Solid dots 
are the standard deviation of Tijfrizi) which is the relative 
deviation of ri from the solid line in panel (a). Note that both 
the spread of z\ and the spread around the solid line vanish as 
1 /%/]V, as if the data were averages of N independent samples. 


t/aexp (((5i - 5j/2)) ^ ya{l - {5j/2)), 


and the ratio of the two different averages becomes 


Vg S 

Va 2yl ’ 


( 21 ) 


which means that the effects discussed here are important 
only when the fluctuations in the data are truly large. 


4- Finite size dependence 

We now examine the spread of Zi and ti , as in 
Fig. 5(b), around the solid line, with special focus on 
how this spread depends on the finite system size. For 
the finite size study we turn to a lower packing fraction, 
(j) = 0.838. The reason for this is that, closer to (/>j (e.g. 
at (/) = 0.840) some configurations for smaller sizes fail to 
reach zero energy in the relaxation step and get jammed 
with z > Ziso, and such events badly complicate the anal¬ 
ysis. 

Fig. 10 which is ti vs zi for cj) = 0.838, the initial shear 
rate 7 = 10“^, and the three sizes, N = 1024, N = 4096, 


and 65536, clearly shows that these data spread more for 
smaller N. Note that the data in Fig. 10 for all different 
sizes have a common behavior, ti pz /r(^i) = Ar(6zi)~^. 
The exponent b = 2.40 is an effective exponent which 
differs from jS/u^ = 2.69 (obtained in Fig. 5) since we 
here make use of data with larger 5zi. 

We introduce three different measures to character¬ 
ize the spread of this data. Two straightforward mea¬ 
sures are s(zi) and s(ti) which are the standard de¬ 
viations of the data. Another measure is the spread 
of Ti away from the line, i.e. the value predicted from 
the known zi, s[ti/ fr{zi)\. These three quantities are 
shown in Fig. 10(b) for number of particles ranging from 
N = 1024 through 65536. To interpret this data we 
first recall that the standard deviation of averages of N 
independent samples is ^ We find that both 

5 ( 2 : 1 ) and s{Ti)/Ta vanish with the exponents —0.54 and 
—0.51 in excellent agreement with this expectation. For 
s[ti/ frizi)] we find a somewhat more complicated be¬ 
haviour with a larger exponent, —0.64, and a question¬ 
able fit to the data. Taken together our data suggest an 
interpretation where both the spread of zi and the spread 
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of Ti around are controlled by independent simple 

stochastic processes. 

IV. DISCUSSION 

The relaxation dynamics around the jamming transi¬ 
tion has been studied before, but with a rather differ¬ 
ent approach [9]: the configurations were first generated 
randomly, then relaxed to a zero-energy state with the 
conjugate gradient method, and after that perturbed by 
a pure affine shear deformation. The relaxation time 
was then determined from the relaxation of such initial 
states by fitting the shear stress to cr((/), t) ~ 
with a = 0.55(5), and the relaxation time was found to 
diverge as r ~ (^j — with ^ = 3.3(1). This expo¬ 
nent is clearly bigger than our /3 « 2.7. One possible 
explanation for this difference is that we in the present 
study get data in the limit of vanishing shear rate in the 
preparation step (i.e. 7 —>■ 0 in the steady state shear¬ 
ing), whereas they in their work apply the pure shear 
deformation suddenly, which is more like a rapid shear¬ 
ing. Indeed, as shown in Fig. 3(a) any given fixed shear 
rate would give too large values for t as one gets close to 
^j, and from analyses of such data one would expect to 
get too high values of the exponent for the divergence. 

We finally want to stress two consequences of the pre¬ 
sented results: We first stress that the above results taken 
together suggest that r is a fundamental quantity that 
controls the overlap J/y and thereby is behind the di¬ 
vergence of other quantites like rjp and rj. For a detailed 
argument we consider the 7 —>■ 0 limit where Tdiss ~ t 
and the N ^ co limit where the spread of zi and ri 
vanish. A given 4> then leads to a well-defined Sz which 
in turn implies a well-defined t and Tdiss ~ t. With the 
additional assumption of a given value for the dimen¬ 
sionless friction, /i = cr/p, power balance between the 
input power Pin = cry ^ /iJy and the dissipated power 
Pdiss = E/Tdiss ~ (5^/Tdiss gives ( 5/7 ~ TdissP- This there¬ 
fore provides a very direct link between the relaxation 
times and rjp ~ S/j. 

Seconly, we note that the relaxation time t we have 
defined here has a different scaling exponent than does 
the time scale associated with rescaling the shear strain 
rate 7 . From Eq. (5) for the 7 —>■ 0 limit and dimensional 
arguments one would expect the deviations due to a finite 
7 to scale as 

~ gUr) ~ gi'i/lScpf), (naive), ( 22 ) 

where the scaling function lima;_>o g{x) = const (for the 
hard disk limit) and the deviations being controlled by 
7 T. This is however not the case. As shown in Eq. (9) the 
data scale with g{'j/\S(j)\^'') where zv = (3 + y, y ~ 1.1, 
which thus is clearly different from the behavior expected 
from dimensional analysis. We hope to be able to return 
to this question elsewhere. 


V. SUMMARY 

To summarize, we have done extensive two-step sim¬ 
ulations, first shearing the system at different constant 
shear rates and then stopping the shearing and letting 
the system relax. At late times of this relaxation, both 
energy and pressure decay exponentially, and we define 
the relaxation time, t, to be the time constant of the ex¬ 
ponential decay of the pressure. We similarly define the 
“dissipation time” from the initial decay immediately af¬ 
ter the shearing is turned off. 

We then show that these two times behave very sim¬ 
ilarly when considering the limit of low shear rates, but 
also that their respective shear rate dependencies are op¬ 
posite. From the expression for Tdiss, Eq- (H), it follows 
immediately that Tdiss diverges with the exponent /3—the 
same divergence as for r]p = p/j —and this is also corrob¬ 
orated by the (/-dependence of t and Tdiss in the small-y 
limit. 

We also show that the relaxation time is directly re¬ 
lated to the lowest vibrational frequency of hard disk 
systems[14], and, furthermore, that this suggests a rela¬ 
tion between t and pp, which should be valid in the small 
7 limit. Fig. 4, provide ample evidence that this actually 
is the case. 

We then turn to a thorough study of the relation be¬ 
tween the contact number and the relaxation time. The 
contact number is a key quantity in the field of jamming 
and we follow Ref. [14] and determine the contact number 
after removing rattlers. With ti and zi from individual 
measurements, ti depends algebraically on the distance 
from isostaticity Szi = Ziso — zi, ti ~ ((5zi)^/“^, with 
/3/u^ « 2.69. 

The same analysis applied to the CDg model gives es¬ 
sentially the same exponent, P/uz ~ 2.63, which pro¬ 
vides additional evidence [17] that the CDq and the RDq 
models are in the same universality class. We consider 
these analysis to be especially robust as they are entirely 
straightforward and do not require data obtained at very 
low shear rates. 

We then turn to effects of the spread of the individual 
Ti for a fixed set of parameters (/, 7 , around its average. 
We first point out that the ordinary arithmetic mean may 
be problematic and that a geometric mean actually in 
some respects works better. We then consider the finite 
size effect where we find that the spread of both the re¬ 
laxation time and the coordination number go as 1 / \/N, 
just as expected for the statistics of N independent vari¬ 
ables. 
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